rm(list = ls())
#working_folder<- "/Users/bgpopescu/Library/CloudStorage/Dropbox/covid_institutions/"
working_folder<-"//Mac/Dropbox-1/covid_paper_replication/"
setwd(working_folder)

#This is to obtain:
#-table_A5

library("readxl")
library("reshape")
library("stringi")
library("stringr")
library("dplyr")
library("tidyr")
library("stargazer")
library("corrplot")
library("Hmisc")
library("ggplot2")
library("glue")
library('gridExtra')
library('grid')
library('lfe')
library('broom')
library("plm")
library("multiwayvcov")
library("lmtest")
library("psych")
library("ggrepel")
library("lemon")
library("sandwich")
library("spdep")
library("xtable")
library('gsynth')
library("pracma")
library("conleyreg")
library("vars")
library("rlist")
library("fastDummies")
library("dplyr")
library("berryFunctions")
library("mediation")
library("sp")
library("pgirmess")
library("fwildclusterboot")
library("geojsonsf")

#Functions

dlshape=function(shploc, shpfile) {
  temp=tempfile()
  download.file(shploc, temp)
  unzip(temp)
  fp <- file.path(temp)
  obj_shp<-rgdal::readOGR(".",shpfile)
  return(obj_shp)}

mean_sd_fun<-function(lm_model) {
  list_empty <-c()
  number_mean<-mean(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  number_sd<-sd(data.frame(lm_model$model)$mean_mobil, na.rm=T)
  list_empty[1]<-number_mean
  list_empty[2]<-number_sd
  names(list_empty)<-c("mean", "sd")
  return(list_empty)
}


week_obs_count<-function(lm_model) {
  list_empty <-c()
  week_no<-length(grep(x = colnames(data.frame(lm_model$model)), pattern = "^week"))
  #Adding 2 because there should be 53 weeks and because I am taking 2 weeks out as in STATA 2 weeks drop
  week_no<-week_no+2
  total_obs<-nobs(lm_model)
  counties<-total_obs/week_no
  list_empty[1]<-week_no
  list_empty[2]<-total_obs
  list_empty[3]<-counties
  names(list_empty)<-c("week_no", "total_obs", "counties")
  return(list_empty)
}


#############################
#Function for calculating SE#
#############################

#Step1. Turning some variables to numeric
data_cov_mob = geojson_sf("./data/data_cov_mob.geojson")
data_cov_mob$pct_manufacturing<-as.numeric(data_cov_mob$pct_manufacturing)

data_cov_mob$kul<-as.numeric(data_cov_mob$kul)
data_cov_mob$nat<-as.numeric(data_cov_mob$nat)
data_cov_mob$frei<-as.numeric(data_cov_mob$frei)
data_cov_mob$soz<-as.numeric(data_cov_mob$soz)
data_cov_mob$pol<-as.numeric(data_cov_mob$pol)
data_cov_mob$inter<-as.numeric(data_cov_mob$inter)


data_cov_mob$gspo<-as.numeric(data_cov_mob$gspo)
data_cov_mob$gkul<-as.numeric(data_cov_mob$gkul)
data_cov_mob$gnat<-as.numeric(data_cov_mob$gnat)

data_cov_mob$gsoz<-as.numeric(data_cov_mob$gsoz)
data_cov_mob$gpol<-as.numeric(data_cov_mob$gpol)
data_cov_mob$ginter<-as.numeric(data_cov_mob$ginter)

#Bridging Networks
data_cov_mob$bridging_total<-data_cov_mob$gspo+ data_cov_mob$gkul+ data_cov_mob$gnat
data_cov_mob$bridging_tot_pop<-data_cov_mob$bridging_total/data_cov_mob$pop_total*1000
summary(data_cov_mob$bridging_tot_pop)
mean_brid<-summary(data_cov_mob$bridging_tot_pop)[4]
data_cov_mob$bridging_tot_pop_dummy<-ifelse(data_cov_mob$bridging_tot_pop>mean_brid, 1, 0)

#Bonding Networks
data_cov_mob$bonding_total<-data_cov_mob$gsoz+ data_cov_mob$gpol+ data_cov_mob$ginter
data_cov_mob$bonding_tot_pop<-data_cov_mob$bonding_total/data_cov_mob$pop_total*1000
summary(data_cov_mob$bonding_tot_pop)
mean_bond<-summary(data_cov_mob$bonding_tot_pop)[4]
data_cov_mob$bonding_tot_pop_dummy<-ifelse(data_cov_mob$bonding_tot_pop>mean_bond, 1, 0)

#All Networks
data_cov_mob$gvereine<-as.numeric(data_cov_mob$gvereine)
data_cov_mob$vereine_tot_pop<-data_cov_mob$gvereine/data_cov_mob$pop_total*1000
summary(data_cov_mob$vereine_tot_pop)
mean_ver<-summary(data_cov_mob$vereine_tot_pop)[4]
data_cov_mob$vereine_tot_pop_dummy<-ifelse(data_cov_mob$vereine_tot_pop>mean_ver, 1, 0)

#Right Wing
data_cov_mob$pct_afd<-as.numeric(data_cov_mob$pct_afd)
summary(data_cov_mob$pct_afd)
mean_pct_afd<-summary(data_cov_mob$pct_afd)[3]
data_cov_mob$pct_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_afd, 1, 0)

#Right Wing Change
data_cov_mob$pct_change_afd<-as.numeric(data_cov_mob$pct_change_afd)
summary(data_cov_mob$pct_change_afd)
mean_pct_change_afd<-summary(data_cov_mob$pct_change_afd)[4]
data_cov_mob$pct_change_afd_dummy<-ifelse(data_cov_mob$pct_afd>mean_pct_change_afd, 1, 0)


table(data_cov_mob$vereine_tot_pop_dummy)
table(data_cov_mob$bridging_tot_pop_dummy)
table(data_cov_mob$pct_afd_dummy)
table(data_cov_mob$pct_change_afd_dummy)


data_cov_mob_week<-subset(data_cov_mob, week=="01")
data_cov_mob_week<-subset(data_cov_mob_week, select = c("GEN", 
                                                        "vereine_tot_pop",
                                                        "bridging_tot_pop", 
                                                        "bonding_tot_pop",
                                                        "pct_afd_dummy",
                                                        "pct_change_afd_dummy"))
summary(data_cov_mob$gender_ratio)
data<-data_cov_mob
data$week_new<-as.numeric(data$week)+1
data$week_new<-paste0("0", data$week_new)  
data$week_new2<-ifelse(nchar(data$week_new)==3, str_remove(data$week_new, "^0+"), data$week_new)
data$week<-data$week_new2
data<-subset(data, select = -c(week_new, week_new2))
data$post_treat<-ifelse(as.numeric(data$week)>10, 1,0)

data$post_treat_vereine_tot_pop_dummy<-data$post_treat*data$vereine_tot_pop_dummy
data$post_treat_bridging_tot_pop_dummy<-data$post_treat*data$bridging_tot_pop_dummy
data$post_treat_bridging_tot_pop<-data$post_treat*data$bridging_tot_pop

data$post_treat_bonding_tot_pop_dummy<-data$post_treat*data$bonding_tot_pop_dummy
data$post_treat_pct_afd_dummy<-data$post_treat*data$pct_afd_dummy
data$post_treat_pct_change_afd_dummy<-data$post_treat*data$pct_change_afd_dummy
data$post_treat_pct_change_afd<-data$post_treat*data$pct_change_afd


#Save to separate file
nuts<-subset(data_cov_mob, week=="00")
nuts2<-subset(nuts, select = c("AGS"))
nuts2$AGS<-as.numeric(nuts2$AGS<-nuts2$AGS)

# Create dummy variables:
dataf <- dummy_cols(data, select_columns = 'AGS')
dataw <- dummy_cols(dataf, select_columns = 'week')
datas <- dummy_cols(dataw, select_columns = 'SN_L')
summary(datas$gender_ratio)


new<-dataw %>% dplyr::select(starts_with("AGS_"))
new2<-dataw %>% dplyr::select(starts_with("week_"))
new3<-datas %>% dplyr::select(starts_with("SN_L"))

list_dummies<-names(new)
#I am removing the counties which get dropped in Stata

#Re-enable this later if necessary
list_dummies<-list_dummies[list_dummies!= "AGS_0"]
list_dummies<-list_dummies[list_dummies!= "AGS_5316"]
list_dummies<-list_dummies[list_dummies!= "AGS_8121"]
list_dummies<-list_dummies[list_dummies!= "AGS_8221"]
list_dummies<-list_dummies[list_dummies!= "AGS_8231"]
list_dummies<-list_dummies[list_dummies!= "AGS_8311"]
list_dummies<-list_dummies[list_dummies!= "AGS_8421"]
list_dummies<-list_dummies[list_dummies!= "AGS_9461"]
list_dummies<-list_dummies[list_dummies!= "AGS_13004"]
list_dummies<-list_dummies[list_dummies!= "AGS_16052"]
list_dummies<-list_dummies[list_dummies!= "AGS_16053"]
list_dummies<-list_dummies[list_dummies!= "AGS_16054"]
list_dummies<-list_dummies[list_dummies!= "AGS_16055"]
list_dummies<-list_dummies[list_dummies!= "AGS_16061"]
list_dummies<-list_dummies[list_dummies!= "AGS_16077"]

list_week_dummies<-names(new2)
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_new2"]
list_week_dummies<-list_week_dummies[list_week_dummies!= "week_numeric"]
list_week_dummies_now1<-list_week_dummies[list_week_dummies!= "week_01"]

#Re-enable this later if necessary
list_week_dummies2<-list_week_dummies[list_week_dummies!= "week_10"]
list_week_dummies2<-list_week_dummies2[list_week_dummies2!= "week_53"]


list_land_dummies<-names(new3)
#Re-enable this later if necessary
list_land_dummies_full<-list_land_dummies[list_land_dummies!= "SN_L"]
list_land_dummies<-list_land_dummies_full[list_land_dummies_full!= "SN_L_01"]
temp_lag<-10


#Creating the interaction terms for bonding
for (i in list_land_dummies){
  for (j in list_week_dummies_now1)
    datas[[paste("inter_", i, "_", j, sep="")]]<-datas[[i]]*datas[[j]]
}

new<-datas %>% dplyr::select(starts_with("inter_"))
list_land_dummies_star<-names(new)

#VEREINE ALL
basic_ver<-c("vereine_tot_pop_dummy", "post_treat", "post_treat_vereine_tot_pop_dummy", "lag_covid_pc")

extended_ver <- c("vereine_tot_pop_dummy", "post_treat", "post_treat_vereine_tot_pop_dummy", "lag_covid_pc",
                  "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                  "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                  "east_ger", "gender_ratio", "pct_students", "uni_pop", "parks_pop")


basic_and_dummies_ver<-list.append(basic_ver, list_dummies, list_week_dummies2)
basic_and_dummies_int_ver<-list.append(basic_ver, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_ver<-list.append(extended_ver, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_ver<-list.append(extended_ver, list_land_dummies, list_week_dummies2, list_land_dummies_star)


spec_basic_ver<-as.formula(paste("mean_mobil~",
                                 paste(basic_and_dummies_ver, collapse= "+")))
spec_basic_int_ver<-as.formula(paste("mean_mobil~",
                                     paste(basic_and_dummies_int_ver, collapse= "+")))

spec_extended_ver<-as.formula(paste("mean_mobil~",
                                    paste(extended_and_week_dummies_ver, collapse= "+")))
spec_extended_int_ver<-as.formula(paste("mean_mobil~",
                                        paste(extended_and_week_dummies_int_ver, collapse= "+")))


#BRIDGING
basic_brid<-c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc")
extended_brid <- c("bridging_tot_pop_dummy", "post_treat", "post_treat_bridging_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students", "uni_pop", "parks_pop")

basic_and_dummies_brid<-list.append(basic_brid, list_dummies, list_week_dummies2)
basic_and_dummies_int_brid<-list.append(basic_brid, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_brid<-list.append(extended_brid, list_land_dummies, list_week_dummies2, list_land_dummies_star)

spec_basic_brid<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_brid, collapse= "+")))
spec_basic_int_brid<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_brid, collapse= "+")))

spec_extended_brid<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_brid, collapse= "+")))
spec_extended_int_brid<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_brid, collapse= "+")))

#BONDING
basic_bond<-c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc")
extended_bond <- c("bonding_tot_pop_dummy", "post_treat", "post_treat_bonding_tot_pop_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                   "east_ger", "gender_ratio", "pct_students", "uni_pop", "parks_pop")

basic_and_dummies_bond<-list.append(basic_bond, list_dummies, list_week_dummies2)
basic_and_dummies_int_bond<-list.append(basic_bond, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_bond<-list.append(extended_bond, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_bond<-list.append(extended_bond, list_land_dummies, list_week_dummies2, list_land_dummies_star)


spec_basic_bond<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_bond, collapse= "+")))
spec_basic_int_bond<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_bond, collapse= "+")))

spec_extended_bond<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_bond, collapse= "+")))
spec_extended_int_bond<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_bond, collapse= "+")))



#PCT. AFD
data$post_treat_pct_afd_dummy
basic_afd<-c("pct_afd_dummy", "post_treat", "post_treat_pct_afd_dummy", "lag_covid_pc")
extended_afd <- c("pct_afd_dummy", "post_treat", "post_treat_pct_afd_dummy", "lag_covid_pc",
                   "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                   "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                  "east_ger", "gender_ratio", "pct_students", "uni_pop", "parks_pop")

basic_and_dummies_afd<-list.append(basic_afd, list_dummies, list_week_dummies2)
basic_and_dummies_int_afd<-list.append(basic_afd, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_afd<-list.append(extended_afd, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_afd<-list.append(extended_afd, list_land_dummies, list_week_dummies2, list_land_dummies_star)

spec_basic_afd<-as.formula(paste("mean_mobil~",
                                  paste(basic_and_dummies_afd, collapse= "+")))
spec_basic_int_afd<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_afd, collapse= "+")))

spec_extended_afd<-as.formula(paste("mean_mobil~",
                                     paste(extended_and_week_dummies_afd, collapse= "+")))
spec_extended_int_afd<-as.formula(paste("mean_mobil~",
                                         paste(extended_and_week_dummies_int_afd, collapse= "+")))




#PCT. AFD CHANGE
data$post_treat_pct_change_afd_dummy
basic_pct_change_afd<-c("pct_change_afd_dummy", "post_treat", "post_treat_pct_change_afd_dummy", "lag_covid_pc")
extended_pct_change_afd <- c("pct_change_afd_dummy", "post_treat", "post_treat_pct_change_afd_dummy", "lag_covid_pc",
                  "pct_turnout", "log_pop_total", "log_gdp_per_cap", "pop_den",
                  "pct_college", "pct_pop_above_60", "pct_pop_under_35", "pct_service_sec_narrow", "pct_manufacturing",
                  "east_ger", "gender_ratio", "pct_students", "uni_pop", "parks_pop")

basic_and_dummies_pct_change_afd<-list.append(basic_pct_change_afd, list_dummies, list_week_dummies2)
basic_and_dummies_int_pct_change_afd<-list.append(basic_pct_change_afd, list_dummies, list_week_dummies2, list_land_dummies_star)

extended_and_week_dummies_pct_change_afd<-list.append(extended_pct_change_afd, list_land_dummies, list_week_dummies2)
extended_and_week_dummies_int_pct_change_afd<-list.append(extended_pct_change_afd, list_land_dummies, list_week_dummies2, list_land_dummies_star)

spec_basic_pct_change_afd<-as.formula(paste("mean_mobil~",
                                 paste(basic_and_dummies_pct_change_afd, collapse= "+")))
spec_basic_int_pct_change_afd<-as.formula(paste("mean_mobil~",
                                     paste(basic_and_dummies_int_pct_change_afd, collapse= "+")))

spec_extended_pct_change_afd<-as.formula(paste("mean_mobil~",
                                    paste(extended_and_week_dummies_pct_change_afd, collapse= "+")))
spec_extended_int_pct_change_afd<-as.formula(paste("mean_mobil~",
                                        paste(extended_and_week_dummies_int_pct_change_afd, collapse= "+")))



temp_lag<-10

#*Model1 - Stata
#*#Very close but slightly different
#xtreg mean_mobil vereine_tot_pop_dummy /*
#  */post_treat /*
#  */post_treat_vereine_tot_pop_dummy /*
#  */lag_covid_pc /*
#  */ags_* /*
#  */week_* /*
#  */i.sn_l##i.week
#Note the variables in Stata which get dropped

datas_city<-subset(datas, BEZ=="Kreisfreie Stadt")
datas_nocity<-subset(datas, BEZ=="Kreis")
mean(datas$pop_den, na.rm=T)
datas_above_meandens<-subset(datas, pop_den>mean(datas$pop_den, na.rm=T))
datas_below_meandens<-subset(datas, pop_den<mean(datas$pop_den, na.rm=T))


mo1_verx<-lm(spec_basic_int_ver, data=datas_city)
summary(mo1_verx)

mo1_ver<-lm(spec_basic_int_ver, data=datas)
summary(mo1_ver)
mo1_ver_df<-broom::tidy(mo1_ver)
mo1_ver_df$rse<-coeftest(mo1_ver, vcov = cluster.vcov(mo1_ver, datas$AGS, force_posdef=T))[, 2]
mo1_ver_df$pv_rse<-coeftest(mo1_ver, vcov = cluster.vcov(mo1_ver, datas$AGS, force_posdef=T))[, 4]
mo1_ver_df2<-mo1_ver_df %>% filter(term %in% basic_ver)


#Model2 - Stata
#Very close but slightly different
#xtreg mean_mobil vereine_tot_pop_dummy /*
#  */post_treat /*
#  */post_treat_vereine_tot_pop_dummy /*
#  */lag_covid_pc /*
#  *ags_*
#  */week_* /*
#  */i.sn_l##i.week /*
# */ pct_turnout log_pop_total log_gdp_per_cap pop_den /*
#  */ pct_college pct_pop_above_60 pct_service_sec_narrow pct_manufacturing

mo8_ver<-lm(spec_extended_int_ver, data=datas)
summary(mo8_ver)
mo8_ver_df<-broom::tidy(mo8_ver)
mo8_ver_df$rse<-coeftest(mo8_ver, vcov = cluster.vcov(mo8_ver, datas$AGS, force_posdef=T))[, 2]
mo8_ver_df$pv_rse<-coeftest(mo8_ver, vcov = cluster.vcov(mo8_ver, datas$AGS, force_posdef=T))[, 4]
mo8_ver_df2<-mo8_ver_df %>% filter(term %in% extended_ver)


#*Model3
#*County FE
#xtreg mean_mobil bridging_tot_pop_dummy /*
#  */post_treat /*
#  */post_treat_bridging_tot_pop_d /*
#  */lag_covid_pc /*
#  */ags_* /*
#  */week_* /*
#  */i.sn_l##i.week


mo1_brid<-lm(spec_basic_int_brid, data=datas)
summary(mo1_brid)
mo1_brid_df<-broom::tidy(mo1_brid)
mo1_brid_df$rse<-coeftest(mo1_brid, vcov = cluster.vcov(mo1_brid, datas$AGS, force_posdef=T))[, 2]
mo1_brid_df$pv_rse<-coeftest(mo1_brid, vcov = cluster.vcov(mo1_brid, datas$AGS, force_posdef=T))[, 4]
mo1_brid_df2<-mo1_brid_df %>% filter(term %in% basic_brid)

basic_and_dummies_int_brid_x<-list.append(basic_brid, list_land_dummies_star)
spec_basic_int_brid_x<-as.formula(paste("mean_mobil~",
                                      paste(basic_and_dummies_int_brid_x, collapse= "+"), "| AGS + week"))

summary(mo1_brid)


# *Model4
# *NO County FE
# *Week FE
# xtreg mean_mobil bridging_tot_pop_dummy /*
#   */post_treat /*
#   */post_treat_bridging_tot_pop_d /*
#   */lag_covid_pc /*
#   *ags_*
#   */week_* /*
#   */i.sn_l##i.week /*
# */ pct_turnout log_pop_total log_gdp_per_cap pop_den /*
#   */ pct_college pct_pop_above_60 pct_service_sec_narrow pct_manufacturing


mo8_brid<-lm(spec_extended_int_brid, data=datas)
summary(mo8_brid)

mo8_brid_df<-broom::tidy(mo8_brid)
mo8_brid_df$rse<-coeftest(mo8_brid, vcov = cluster.vcov(mo8_brid, datas$AGS, force_posdef=T))[, 2]
mo8_brid_df$pv_rse<-coeftest(mo8_brid, vcov = cluster.vcov(mo8_brid, datas$AGS, force_posdef=T))[, 4]
mo8_brid_df2<-mo8_brid_df %>% filter(term %in% extended_brid)


# *Model5
# *County FE
# xtreg mean_mobil bonding_tot_pop_dummy /*
#   */post_treat /*
#   */post_treat_bonding_tot_pop_d /*
#   */lag_covid_pc /*
#   */ags_* /*
#   */week_* /*
#   */i.sn_l##i.week


mo1_bond<-lm(spec_basic_int_bond, data=datas)
summary(mo1_bond)
mo1_bond_df<-broom::tidy(mo1_bond)
mo1_bond_df$rse<-coeftest(mo1_bond, vcov = cluster.vcov(mo1_bond, datas$AGS, force_posdef=T))[, 2]
mo1_bond_df$pv_rse<-coeftest(mo1_bond, vcov = cluster.vcov(mo1_bond, datas$AGS, force_posdef=T))[, 4]
mo1_bond_df2<-mo1_bond_df %>% filter(term %in% basic_bond)


# *Model6
# *NO County FE
# *Week FE
# xtreg mean_mobil bonding_tot_pop_dummy /*
#   */post_treat /*
#   */post_treat_bonding_tot_pop_d /*
#   */lag_covid_pc /*
#   *ags_*
#   */week_* /*
#   */i.sn_l##i.week /*
# */ pct_turnout log_pop_total log_gdp_per_cap pop_den /*
#   */ pct_college pct_pop_above_60 pct_service_sec_narrow pct_manufacturing

mo8_bond<-lm(spec_extended_int_bond, data=datas)
summary(mo8_bond)
mo8_bond_df<-broom::tidy(mo8_bond)
mo8_bond_df$rse<-coeftest(mo8_bond, vcov = cluster.vcov(mo8_bond, datas$AGS, force_posdef=T))[, 2]
mo8_bond_df$pv_rse<-coeftest(mo8_bond, vcov = cluster.vcov(mo8_bond, datas$AGS, force_posdef=T))[, 4]
mo8_bond_df2<-mo8_bond_df %>% filter(term %in% extended_bond)



model1<-mo1_ver_df2
model2<-mo8_ver_df2
model3<-mo1_brid_df2
model4<-mo8_brid_df2
model5<-mo1_bond_df2
model6<-mo8_bond_df2


model1$order<-NA
model2$order<-NA
model3$order<-NA
model4$order<-NA
model5$order<-NA
model6$order<-NA

model1$order[model1$term=="vereine_tot_pop_dummy"]<-1
model2$order[model2$term=="vereine_tot_pop_dummy"]<-1
model3$order[model3$term=="bridging_tot_pop_dummy"]<-2
model4$order[model4$term=="bridging_tot_pop_dummy"]<-2
model5$order[model5$term=="bonding_tot_pop_dummy"]<-3
model6$order[model6$term=="bonding_tot_pop_dummy"]<-3

model1$order[model1$term=="post_treat"]<-4
model2$order[model2$term=="post_treat"]<-4
model3$order[model3$term=="post_treat"]<-4
model4$order[model4$term=="post_treat"]<-4
model5$order[model5$term=="post_treat"]<-4
model6$order[model6$term=="post_treat"]<-4

model1$order[model1$term=="post_treat_vereine_tot_pop_dummy"]<-5
model2$order[model2$term=="post_treat_vereine_tot_pop_dummy"]<-5
model3$order[model3$term=="post_treat_bridging_tot_pop_dummy"]<-6
model4$order[model4$term=="post_treat_bridging_tot_pop_dummy"]<-6
model5$order[model5$term=="post_treat_bonding_tot_pop_dummy"]<-7
model6$order[model6$term=="post_treat_bonding_tot_pop_dummy"]<-7


model2$order[model2$term=="pct_turnout"]<-8
model4$order[model4$term=="pct_turnout"]<-8
model6$order[model6$term=="pct_turnout"]<-8

model2$order[model2$term=="log_pop_total"]<-9
model4$order[model4$term=="log_pop_total"]<-9
model6$order[model6$term=="log_pop_total"]<-9

model2$order[model2$term=="log_gdp_per_cap"]<-10
model4$order[model4$term=="log_gdp_per_cap"]<-10
model6$order[model6$term=="log_gdp_per_cap"]<-10

model2$order[model2$term=="pop_den"]<-11
model4$order[model4$term=="pop_den"]<-11
model6$order[model6$term=="pop_den"]<-11

model2$order[model2$term=="pct_college"]<-12
model4$order[model4$term=="pct_college"]<-12
model6$order[model6$term=="pct_college"]<-12

model2$order[model2$term=="pct_pop_above_60"]<-13
model4$order[model4$term=="pct_pop_above_60"]<-13
model6$order[model6$term=="pct_pop_above_60"]<-13

model2$order[model2$term=="pct_pop_under_35"]<-14
model4$order[model4$term=="pct_pop_under_35"]<-14
model6$order[model6$term=="pct_pop_under_35"]<-14

model2$order[model2$term=="pct_service_sec_narrow"]<-15
model4$order[model4$term=="pct_service_sec_narrow"]<-15
model6$order[model6$term=="pct_service_sec_narrow"]<-15

model2$order[model2$term=="pct_manufacturing"]<-16
model4$order[model4$term=="pct_manufacturing"]<-16
model6$order[model6$term=="pct_manufacturing"]<-16

model2$order[model2$term=="east_ger"]<-17
model4$order[model4$term=="east_ger"]<-17
model6$order[model6$term=="east_ger"]<-17

model1$order[model1$term=="lag_covid_pc"]<-18
model2$order[model2$term=="lag_covid_pc"]<-18
model3$order[model3$term=="lag_covid_pc"]<-18
model4$order[model4$term=="lag_covid_pc"]<-18
model5$order[model5$term=="lag_covid_pc"]<-18
model6$order[model6$term=="lag_covid_pc"]<-18


model2$order[model2$term=="gender_ratio"]<-19
model4$order[model4$term=="gender_ratio"]<-19
model6$order[model6$term=="gender_ratio"]<-19

model2$order[model2$term=="pct_students"]<-20
model4$order[model4$term=="pct_students"]<-20
model6$order[model6$term=="pct_students"]<-20

model2$order[model2$term=="uni_pop"]<-21
model4$order[model4$term=="uni_pop"]<-21
model6$order[model6$term=="uni_pop"]<-21

model2$order[model2$term=="parks_pop"]<-22
model4$order[model4$term=="parks_pop"]<-22
model6$order[model6$term=="parks_pop"]<-22


model1$estimate<-round(as.numeric(model1$estimate), 3)
model2$estimate<-round(as.numeric(model2$estimate), 3)
model3$estimate<-round(as.numeric(model3$estimate), 3)
model4$estimate<-round(as.numeric(model4$estimate), 3)
model5$estimate<-round(as.numeric(model5$estimate), 3)
model6$estimate<-round(as.numeric(model6$estimate), 3)


model1$std_er_star<-paste("(", round(as.numeric(model1$rse),3), ")",
                          ifelse(as.numeric(model1$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model1$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model1$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model2$std_er_star<-paste("(", round(as.numeric(model2$rse),3), ")",
                          ifelse(as.numeric(model2$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model2$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model2$pv_rse)<0.05,"*",
                                               ""))), sep = "")
model3$std_er_star<-paste("(", round(as.numeric(model3$rse),3), ")",
                          ifelse(as.numeric(model3$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model3$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model3$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model4$std_er_star<-paste("(", round(as.numeric(model4$rse),3), ")",
                          ifelse(as.numeric(model4$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model4$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model4$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model5$std_er_star<-paste("(", round(as.numeric(model5$rse),3), ")",
                          ifelse(as.numeric(model5$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model5$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model5$pv_rse)<0.05,"*",
                                               ""))), sep = "")

model6$std_er_star<-paste("(", round(as.numeric(model6$rse),3), ")",
                          ifelse(as.numeric(model6$pv_rse)<0.001,"***",
                                 ifelse(as.numeric(model6$pv_rse<0.01),"**",
                                        ifelse(as.numeric(model6$pv_rse)<0.05,"*",
                                               ""))), sep = "")


mode1_b<-subset(model1, select = c("order", "term", "estimate", "std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_ver)[1], 3), "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_ver)[2], 3), "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_ver)[3], 3), "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_ver)[1], 3), "", "")
county_fe<-c(27, "county_fe", "Yes", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_ver)$adj.r.squared, 3), "", "")
obs<-c(31, "obs", nobs(mo1_ver), "", "")

mode1_b<-rbind(mode1_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode1_b$model<-"Model1"
mode1_b$order<-as.numeric(mode1_b$order)
names(mode1_b)<-paste0(names(mode1_b), "_1")


mode2_b<-subset(model2, select = c("order", "term", "estimate", "std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_ver)[1], 3), "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_ver)[2], 3), "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_ver)[3], 3), "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_ver)[1], 3), "", "")
county_fe<-c(27, "county_fe", "No", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_ver)$adj.r.squared, 3), "", "")
obs<-c(31, "obs", nobs(mo8_ver), "", "")

mode2_b<-rbind(mode2_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode2_b$model<-"Model2"
mode2_b$order<-as.numeric(mode2_b$order)
names(mode2_b)<-paste0(names(mode2_b), "_2")


mode3_b<-subset(model3, select = c("order", "term", "estimate", "std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_brid)[1], 3), "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_brid)[2], 3), "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_brid)[3], 3), "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_brid)[1], 3), "", "")
county_fe<-c(27, "county_fe", "Yes", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_brid)$adj.r.squared, 3), "", "")
obs<-c(31, "obs", nobs(mo1_brid), "", "")

mode3_b<-rbind(mode3_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode3_b$model<-"Model3"
mode3_b$order<-as.numeric(mode3_b$order)
names(mode3_b)<-paste0(names(mode3_b), "_3")


mode4_b<-subset(model4, select = c("order", "term", "estimate", "std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_brid)[1], 3), "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_brid)[2], 3), "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_brid)[3], 3), "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_brid)[1], 3), "", "")
county_fe<-c(27, "county_fe", "No", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_brid)$adj.r.squared, 3), "", "")
obs<-c(31, "obs", nobs(mo8_brid), "", "")

mode4_b<-rbind(mode4_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode4_b$model<-"Model4"
mode4_b$order<-as.numeric(mode4_b$order)
names(mode4_b)<-paste0(names(mode4_b), "_4")


mode5_b<-subset(model5, select = c("order", "term", "estimate", "std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo1_bond)[1], 3), "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo1_bond)[2], 3), "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo1_bond)[3], 3), "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo1_bond)[1], 3), "", "")
county_fe<-c(27, "county_fe", "Yes", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo1_bond)$adj.r.squared, 3), "", "")
obs<-c(31, "obs", nobs(mo1_bond), "", "")

mode5_b<-rbind(mode5_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode5_b$model<-"Model5"
mode5_b$order<-as.numeric(mode5_b$order)
names(mode5_b)<-paste0(names(mode5_b), "_5")

mode6_b<-subset(model6, select = c("order", "term", "estimate", "std_er_star"))
means<-c(23, "mean_mob", round(mean_sd_fun(mo8_bond)[1], 3), "", "")
sds<-c(24, "sd_mob", round(mean_sd_fun(mo8_bond)[2], 3), "", "")
no_cross_sections<-c(25, "cross_sections", round(week_obs_count(mo8_bond)[3], 3), "", "")
no_time_units<-c(26, "time_units", round(week_obs_count(mo8_bond)[1], 3), "", "")
county_fe<-c(27, "county_fe", "No", "", "")
week_fe<-c(28, "week_fe", "Yes", "", "")
weekxland_fe<-c(29, "weekxland_fe", "Yes", "", "")
adj_rsq<-c(30, "adj_rsq", round(summary(mo8_bond)$adj.r.squared, 3), "", "")
obs<-c(31, "obs", nobs(mo8_bond), "", "")

mode6_b<-rbind(mode6_b, means, sds, no_cross_sections, no_time_units, county_fe,
               week_fe, weekxland_fe, adj_rsq, obs)
mode6_b$model<-"Model6"
mode6_b$order<-as.numeric(mode6_b$order)
names(mode6_b)<-paste0(names(mode6_b), "_6")


#Step1: Collecting the models
list_order<-c(mode1_b$order_1, mode2_b$order_2, mode3_b$order_3,
              mode4_b$order_4, mode5_b$order_5, mode6_b$order_6)
#Step2: Ordering rows
list_order<-unique(sort(list_order))
list_order_df<-data.frame(list_order)
names(list_order_df)<-"order"
list_order_df$order<-as.numeric(list_order_df$order)


new_df2<-left_join(list_order_df, as.data.frame(mode1_b), by = c("order"="order_1"))


new_df3<-left_join(new_df2, as.data.frame(mode2_b), by = c("order"="order_2"))
new_df3$model<-new_df3$model_1
new_df3$model<-ifelse(is.na(new_df3$model), new_df3$model_2, new_df3$model)
new_df3$term<-new_df3$term_1
new_df3$term<-ifelse(is.na(new_df3$term_1), new_df3$term_2, new_df3$term)
new_df3<-subset(new_df3, select = -c(term_1, model_1, term_2, model_2))

new_df4<-left_join(new_df3, as.data.frame(mode3_b), by = c("order"="order_3"))
new_df4$model<-ifelse(is.na(new_df4$model), new_df4$model_3, new_df4$model)
new_df4$term<-ifelse(is.na(new_df4$term), new_df4$term_3, new_df4$term)
new_df4<-subset(new_df4, select = -c(term_3, model_3))

new_df5<-left_join(new_df4, as.data.frame(mode4_b), by = c("order"="order_4"))
new_df5$model<-ifelse(is.na(new_df5$model), new_df5$model_4, new_df5$model)
new_df5$term<-ifelse(is.na(new_df5$term), new_df5$term_4, new_df5$term)
new_df5<-subset(new_df5, select = -c(term_4, model_4))

new_df6<-left_join(new_df5, as.data.frame(mode5_b), by = c("order"="order_5"))
new_df6$model<-ifelse(is.na(new_df6$model), new_df6$model_5, new_df6$model)
new_df6$term<-ifelse(is.na(new_df6$term), new_df6$term_5, new_df6$term)
new_df6<-subset(new_df6, select = -c(term_5, model_5))

new_df7<-left_join(new_df6, as.data.frame(mode6_b), by = c("order"="order_6"))
new_df7$model<-ifelse(is.na(new_df7$model), new_df7$model_6, new_df7$model)
new_df7$term<-ifelse(is.na(new_df7$term), new_df7$term_6, new_df7$term)
new_df7<-subset(new_df7, select = -c(term_6, model_6))


new_df7$formal<-NA
new_df7$formal[new_df7$term=="vereine_tot_pop_dummy"]<-"All Networks"
new_df7$formal[new_df7$term=="bridging_tot_pop_dummy"]<-"Bridging"
new_df7$formal[new_df7$term=="bonding_tot_pop_dummy"]<-"Bonding"
new_df7$formal[new_df7$term=="post_treat"]<-"Post Week 11 Dummy"
new_df7$formal[new_df7$term=="post_treat_vereine_tot_pop_dummy"]<-"Post Week 11 x All Networks"
new_df7$formal[new_df7$term=="post_treat_bridging_tot_pop_dummy"]<-"Post Week 11 x Bridging"
new_df7$formal[new_df7$term=="post_treat_bonding_tot_pop_dummy"]<-"Post Week 11 x Bonding"
new_df7$formal[new_df7$term=="lag_covid_pc"]<-"Lag Covid per Cap."
new_df7$formal[new_df7$term=="pct_turnout"]<-"Pct. Turnout"
new_df7$formal[new_df7$term=="log_pop_total"]<-"Log Pop. Total"
new_df7$formal[new_df7$term=="log_gdp_per_cap"]<-"Log GDP/capita"
new_df7$formal[new_df7$term=="pop_den"]<-"Pop. Density"
new_df7$formal[new_df7$term=="pct_college"]<-"Pct. College"
new_df7$formal[new_df7$term=="pct_pop_above_60"]<-"Pct. Pop. above 60"
new_df7$formal[new_df7$term=="pct_pop_under_35"]<-"Pct. Pop. below 35"
new_df7$formal[new_df7$term=="pct_service_sec_narrow"]<-"Pct. Empl. Hospitality and Transport"
new_df7$formal[new_df7$term=="pct_manufacturing"]<-"Pct. Empl. Manufacturing"
new_df7$formal[new_df7$term=="east_ger"]<-"East Germany"
new_df7$formal[new_df7$term=="gender_ratio"]<-"Gender Ratio"
new_df7$formal[new_df7$term=="pct_students"]<-"Pct. Students"
new_df7$formal[new_df7$term=="uni_pop"]<-"University-Population Ratio"
new_df7$formal[new_df7$term=="parks_pop"]<-"Parks-Population Ratio"
new_df7$formal[new_df7$term=="mean_mob"]<-"Mean Mobility"
new_df7$formal[new_df7$term=="sd_mob"]<-"SD Mobility"
new_df7$formal[new_df7$term=="cross_sections"]<-"No. Cross-Sections"
new_df7$formal[new_df7$term=="time_units"]<-"No. Time Units"
new_df7$formal[new_df7$term=="county_fe"]<-"County FE"
new_df7$formal[new_df7$term=="week_fe"]<-"Week FE"
new_df7$formal[new_df7$term=="weekxland_fe"]<-"Week X Land FE"
new_df7$formal[new_df7$term=="adj_rsq"]<-"Adj. R sq."
new_df7$formal[new_df7$term=="obs"]<-"Observations"



new_df8<-subset(new_df7, select = c(formal, estimate_1, std_er_star_1,
                                    estimate_2, std_er_star_2,
                                    estimate_3, std_er_star_3,
                                    estimate_4, std_er_star_4,
                                    estimate_5, std_er_star_5,
                                    estimate_6, std_er_star_6))


#Step3 Order columns
#final1<-new_df8[ , order(names(new_df8))]
final1<-new_df8
new_row<-c("Variable", rep(c("Estimate", "Robust SE"), 6))
final2<-rbind(new_row, final1)


#Step7 centering coefficients
object_latex<-xtable(final2, type = "latex", 
                     caption = "Results",
                     #Note there are X columns to be aligned below.
                     #This is because of the index column.
                     align = "llllllllllllll")


comment <- list(pos = list(0), command = NULL)
comment$pos[[1]] <- c(nrow(object_latex))
comment$command <- c(paste("\\bottomrule\n",
                           "\\multicolumn{13}{l} {\\parbox[t]{35cm}\n",
                           "{\\footnotesize \\textit{Note}: 
This table shows the effect of networks on mobility in the post-lockdown period. 
The dependent variable is mobility as measured by the government. 
The interactions between the period after week 11 and different types of networks (all networks, bridging, bonding)
are the variables of interest. All networks, bridging or bonding take the value of 1 for counties with 
density of these individual networks above the mean. Similarly, the post-lockdown period is marked with 1.
The standard errors are clustered at a county level in round brackets. Significant at ∗10\\%, ∗∗5\\%, and ∗∗∗1\\%.}} \\\\ \n", sep = ""))



#Step8 Printing latex tex
tableLines<-print(object_latex,
                  #The following line controls where the lines are drawn
                  hline.after=c(-1, 1, 23),
                  floating = FALSE,
                  file = "./paper/tables/table_A5.tex", 
                  #caption.placement = "top",
                  add.to.row = comment,
                  booktabs = TRUE,
                  include.colnames=F,
                  include.rownames=FALSE,
                  sanitize.text.function=function(x){x})

multicolumns <- "& \\\\multicolumn{2}{c}{Model 1} & \\\\multicolumn{2}{c}{Model 2} & 
\\\\multicolumn{2}{c}{Model 3} & \\\\multicolumn{2}{c}{Model 4} & \\\\multicolumn{2}{c}{Model 5}
& \\\\multicolumn{2}{c}{Model 6}\\\\\\\\"
clines<-"\\\\cmidrule(lr){2-3} \\\\cmidrule(lr){4-5} \\\\cmidrule(lr){6-7} \\\\cmidrule(lr){8-9}
\\\\cmidrule(lr){10-11} \\\\cmidrule(lr){12-13}\\\\\\\\"
tableLines <- sub ("\\\\toprule\\n", paste0 ("\\\\toprule\n", multicolumns, clines, "\n"), tableLines) ## booktabs = TRUE
#tableLines <- sub ("\\\\hline\\n",   paste0 ("\\\\hline\n",   multicolumns, "\n"), tableLines) ## booktabs = FALSE
writeLines(tableLines, con = "./paper/tables/table_A5.tex")
